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Abstract. We describe a method for determining the 
limb polarization and limb darkening of stars in eclipsing 
binary systems, by inverting photometric and polarimetric 
light curves. 

Because of the ill-conditioning of the problem, we use 
the Backus-Gilbert method to control the resolution and 
stability of the recovered solution, and to make quantita- 
tive estimates of the maximum accuracy possible. Using 
this method we confirm that the limb polarization can in- 
deed be recovered, and demonstrate this with simulated 
data, thus determining the level of observational accuracy 
required to achieve a given accuracy of reconstruction. 
This allows us to set out an optimal observational strat- 
egy, and to critcally assess the claimed detection of limb 
polarization in the Algol system. 

The use of polarization in stars has been proposed as a 
diagnostic tool in microlensing surveys by Simmons et al. 
( 1995 ), and we discuss the extension of this work to the 
case of microlensing of extended sources. 

Key words: Polarization - Methods: data analysis - 
Stars: binaries: eclipsing - Cosmology: gravitational tens- 
ing 



in principle be detectable. Theoretical polarimetric light 
curves for this situation have been calculated by Landi 



1. Introduction 

Scattered light emerging from a stellar atmosphere is ex- 
pected to be partially linearly polarized. This effect should 
be greatest at the limb of the star, with a pure electron 
scattering atmosphere giving a limb polarization of 11.7% 
(Chandrasekhar 1946, 195C). Less idealised calculations 



(Collins & Buerghcr 1974) suggest a lower degree of limb 
polarization, around 2% for early-type stars, though this 
figure is rather sensitive to the ionization state of the out- 
ermost atmospheric layers. Clearly a spherically symmet- 
ric star will exhibit no net polarization, but if the sym- 
metry is broken by an eclipse the limb polarization should 



DegPInnocenti et al. (g_98S|) . 

The first detection of this 'Chandrasekhar effect' was 



Send offprint requests to: N Gray 

* Present address: British Antarctic Survey, Cambridge, CB3 
OET, UK 



in the Algol system (Kemp et al. 1983). This data was 
analysed by Wilson and Liou (1993), but the complex- 
ity of the system and the amount of modelling involved 
in the analysis prevented them from making any reliable 
estimate of the limb polarization. 

The inversion of polarimetric light curves from eclips- 
ing binary stars should allow the limb polarization of the 
eclipsed star to be measured (and the photometric light 
curve should similarly give the limb darkening). In fact 
this inverse problem is highly ill-conditioned, and relat- 
ing observations to stellar atmosphere models is therefore 
far from straightforward. In this paper we investigate the 
practical feasibility of determining of limb polarization by 
this method. 

This allows us to address three closely related issues: 

1. We develop a method of obtaining the polarization at 
a point on the stellar disc, and of estimating the er- 
ror on this value. This is based on the Backus-Gilbert 
inversion technique. 

2. We determine the maximium accuracy possible in de- 
termining limb polarization, given a number of data 
points and a noise level. 

3. Thus, we are able to put forward an observational 
strategy which should allow the best measurement of 
limb polarization. 

In section 2 we give a brief overview of the problem. 
We set out the formalism of the Backus-Gilbert method 
in section 3, and discuss its suitability for the problem at 
hand. Section 4 contains the calculations for the specific 
case of eclipsing binary stars, and section 5 presents the 
results of the inversion scheme when applied to simulated 
data, and the conclusions that can be drawn from these. 
Section 6 considers a simplified analogue of the Algol sys- 
tem, comparing the theoretical polarization profile with 
the best resolution current measurements can achieve. 
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2. An overview of the problem 

The radiation field across the stellar disk can be described 
by the stokes parameters, which give the unpolarized and 
polarized intensity. These are usually denoted by /, Q, 
U, and V (these p arame ters are discussed more fully in 
Clarke & Grainger ( 1971 )). Here we shall take the circular 
polarization V to be zero. U and Q give the state of linear 
polarization of the radiation. Q is obtained by measuring 
the difference in intensity in two perpendicular directions, 
and U by measuring the difference when the polarimiter 
is rotated through 45°. The degree of linear polarization 

1 /2 

may simply be written as (Q 2 + U 2 ) /I and the position 
angle as A arctan U/Q. Under rotation of axes by t/j, Q and 
U simply transform as t/ ncw = Cold cos 2'0 — Qo\d sm 20 
and Qnow = U \a sin 2tp + Q \a cos 2ip. I is invariant (as is 
V). ' 

We will assume that the light emerging from the stellar 
disk is partially linearly polarized perpendicular to the 
radial direction. One can calculate the total flux from the 
eclipsed star by integrating the intensity / over the visible 
part of the disc. Similarly, to obtain the total polarized 
flux one simply integrates the Stokes parameter Q over 
the same part of the disc, introducing a rotation factor to 
transform the polarized intensity at each point into the 
appropriate reference frame. 



I(r)A(r,cf),t)rdrd(j> (1) 



disc 



Q(r) A(r, <f>,t) cos 20r drd0 (2) 



disc 



1 

R 2 



Q(r) A(r,(f>,t) sin20rdrd</> (3) 



where (r, (f>) are polar coordinates on the stellar disk, R is 
the distance from the observer to the star, and A(r, <j>, t) 
is zero for an occulted point on the disk, 1 otherwise. 

In this paper we address the problem of extracting 
I(r) and Q(r) from measurements of F\, Fq and F\j dur- 
ing an eclipse. One might attempt to do this by modelling 
I(r) and Q(r) for the stellar atmosphere, performing the 
integrals above and using the data to fit the free param- 
eters in the stellar atmosphere model. But the problem 
must be treated more carefully. It is a general property 
of inverse problems, such as the present case, that a very 
large range of source models are consistent with a given 
data set. This is essentially due to the smoothing prop- 
erties of the kernel (here, the occultation function A). 
Consequently, a forward-modelling approach can give ap- 
parently accurate, but actually highly misleading, results. 
A thorough treatment of inverse problems in astronomy 
can be found in Craig & Brown (1986); here we will sim- 



ply state that meaningful results can only be obtained by 
controlling the instabilities caused by discrete, noisy data. 
We have used the Backus-Gilbert inversion technique to 



achieve this control. The principal advantage of the Ba- 
ckus-Gilbert method from our point of view is that it al- 
lows us a qualitative understanding of, and thus an ex- 
plicit quantitative control over, the compromise between 
bias and stability in our inversion. 

3. The Backus-Gilbert method 



The Backus-Gilbert method is one of a family of methods 
for attacking inverse problems. Ther e are several general 
introductions to the method (Park er (|1977 ) gives an excel- 
lent review, and Loredo & Epstein ( 198S ) gives an example 
in an astrophysical context) - we give a minimal introduc- 
tion here, partly in order to fix our notation. There is a 
more formal discussion of the method in an appendix, and 
a thor ough treatment can be found in Backus & Gilbert 
( |1970D . 

We wish to recover an underlying function u(r), in this 
case, either an intensity, or the run of polarization with 
radius; the parameter r here is the projected radial dis- 
tance from the centre of the eclipsed star's disk, in units 
where the star's radius is r = 1. We cannot measure this 
underlying function directly, but instead measure an inte- 
gral of it, F(s), which might be the observed luminosity 
or polarisation of the eclipsed star. This is related to the 
underlying function through 



F(s) = [ u(r)K(r;s)dr, 
Jo 



(4) 



where s is the vector between the centres of the star and 
its occultor, and the kernel K(r; s) can be calculated a 
priori. The set of observations that we make, fi = F(si), 
is therefore related to this underlying function by 



fi 



u(r)Ki(r) dr + »j 



(5) 



where Ki(r) = K(r; Si) and Ui is a random admixture of 
noise. From these measurements we wish to produce an 
estimator u(r) of the underlying function u(r). 

For the Backus-Gilbert method, we suppose that the 
mean of the estimator and the underlying function are 
related by an averaging kernel A(r, r'), through 



E{u{r)) 



A(r, r')u(r') dr' . 



(G) 



Since we do not know the underlying function, the aver- 
aging kernel is of no use to us directly; however we can 
study its properties, and use our data fi in such a way 
as to optimise those properties, and so minimise the de- 
pendence of the estimate u{r) on the underlying function 
and the noise. It is clear from Eqn. ^) that u(r) = u{r) 
if A(r, r') is the Dirac delta function. However, such a so- 
lution is highly sensitive to noise in the data and hence 
very unstable. We attain stability by increasing the width 
of A(r, r'), and hence smoothing the recovered value over 
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a wider region of the source - that is, we minimise the 
width of A(r, r') subject to the conflicting demand that 
the resolution of A(r, r') is sufficient for the problem at 
hand, in a sense we shall make more precise below. 

We relate our data /j to our estimator u(r) through a 
set of response kernels qi(r), which produce an estimate 
of the underlying function through 



(7) 



If we substitute Eqn. (|^) into Eqn. (^), and assume 
E QDj Qi{ r ) n i) — 0> then we can compare with Eqn. (0), 
and obtain 



A(r,r') = J2^(r)K t (r'). 



(8) 



This allows us to form some measure of the width of 
A(r, r') such as 



A = j (r-r')'[A(r,r')] 2 dr', 



(9) 



qz(r)W tj (r)q 3 (r) 



= q(r) T W(r)q(r), 

which depends on qi, and depends on K% through the def- 
inition 



W l3 {r)= / (r' -r) 2 K t (r')K 3 {r')dr' 



(10) 



This is the standard definition of the width; others are rea- 
sonable, and may be preferable in different circumstances. 

Again assuming that E (q ■ n) =0, we can also form a 
measure of the stability of Eqn. (Q) 

B = Varw(r) = q(r) T Sq(r), (11) 

which depends on qi and the noise covariance matrix SVj = 
E(niHj). In our simulations below, we take the Hi to be 
independent and Gaussian with standard deviation a; in 
this case, Su = 5ijO~ 2 , and our assumption E (q ■ n) = 
is true. 

Finally, the demand that A(r, r') in Eqn. (^) have unit 
area, leads to the constraint 



q(r)-R=l, 



(12) 



where Ri = J Ki(r) dr. 

The Backus-Gilbert method consists of finding those 
qi(r) which minimise 

A + XB = q(r) ■ [W(r) + AS] • q(r) 

(r-r') 2 [A(r, r')} 2 dr' + A Var u(r), 



for some selected parameter A, subject to the constraint 
q ■ R = 1 . The minimisation problem has explicit analytic 
solutions 



9a0) 



[Wj^ + XS]- 1 R 
R ■ [W(r) + AS*] -1 ■ R 



(13) 



in terms of the parameter A, and these different solutions, 
when combined with the data ft using Eqn. (Q), give dif- 
ferent estimators u\(r). The nature of the trade-off in 
the minimisation is clear: in order to improve the stabil- 
ity of the recovery, we choose a A which makes A(r, r') 
broader, and so generate response kernels qi which extend 
the weighted average over a greater number of the data 
points /,-. The cost of this is that the estimate of the re- 
covered point will be biased by the inclusion of the extra 
data, and this will be more marked when the underlying 
function is rapidly varying. 

The first important point about the Backus-Gilbert 
method is that the parameter A allows us to adjudicate 
between the conflicting demands of minimising the width 
of the kernel A(r, r') and minimising the sensitivity of the 
recovered value (which is a realisation of the statistical 
variable u\(r)) to the measurement noise, and that this 
adjudication can be done prior to any data being collected, 
based only on the characteristics of the kernel K(r;s) and 
the noise. 

Secondly, we must emphasise that the q x (r) we obtain 
gives us, through Eqn. (Q), a single point in the recovered 
function, u\ (r) . This means that in this simplest version of 
the Backus-Gilbert method we must perform the inversion 
for each value of r for which we wish to find u\(r). Since 
the calculation of the coefficients q A (r) involves a matrix 
inversion, which is an n 3 procedure, it can be computa- 
tionally expensive, but this limitation is acceptable in our 
particular case, as we only wish to recover the polarisa- 
tions at the limb, r = 1. This feature has the compen- 
sation that we can if necessary select a different optimal 
value of A for each recovered point. 

4. Kernels for the eclipsing case 

In this section we derive the kernels for the case when 
one star is parti ally eclipsed b y another (a future paper 



(Coleman et al. frn preparation ) will discuss the case of a 
gravitational lens). The geometry of this situation is as 
shown in Fig. [|. We can take the occultor to be opaque, 
so that its 'transfer function' Ao(ro), where ro is the (pro- 
jected) radius from the centre of the occultor, is 



A (r ) 



for ?'o < p 

1 for r > p 



(14) 



This makes A(r, x) s), the kernel in terms of r, a function 
of ip(r;s). In terms of i(r), the intensity of the star as 
a function of (projected) radius, the total flux from the 
eclipsed star is 



(15) 



Fj(s) = / i(r)A(r, x, s) (d 2 r = rdrdx). 

J area 

Thus define A(r; s) — r A(r, %; s) d%, so that 



Fi(s) = / i(r)A(r;s)dr, 
Jo 



(16) 
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we therefore find that the total polarized flux in the Q di- 
rection, measured when the centres are a distance s apart, 
is 

F Q (s(t),<t>(t))= [ P(r)i Q (r; S ,0)dr, (20) 
Jo 

and similarly for A\j (r, \\ s , 4>) and F\j(s,4>). This is now 
in the form of Eqn. ([|). Setting I — s cos cf>, we can thus 
see that 



-4q0;s,< 



Au(r;s, ( 



r cos 2cf)(t) sin ip(r; s(t)) 



I 2 



2r ( 2- - 1 ) ly /l 



r sin 2<j) sin if) 



I 



2 \ 2 



= 4r- 1 - -y 7V 7 ! - 



r 



(21) 



(22) 



These kernels are broad and smooth, hence the ill- 
conditioning of the inverse problem. 



Fig. 1. The eclipsed star has radius 1, and the occultor ra- 
dius p. (r, x) are polar coordinates in the projected plane 
of the source, and (s(t), (p(t)) are the coordinates of the 
centre of the occultor; both \ and cf> are taken from the 
radius which touches the point of closest approach. When 
the centres of the two stars are a distance s apart, the oc- 
cultor cuts off an angle tp(r; s, p) of an annulus of radius r, 
centred on the eclipsed star. 



putting Eqn. (|1 

It is easy to see that A(r;s) 
7 = cos[^(r; s)/2] we find that 

^ 7 = (r 2 + s 2 - p 2 )/2rs 



into the form of Eqn. 

2-7T — 1p~{s) 



'Jo 



1 d%- Writing 



7 



-1 



1 < 7 
-1 < 7 < 1 
7 < -1 



(17) 



A(r;s) 



P 



(18) 



We can now write 

2r(n - tp/2) 

= 2rarccos(— 7) , \s — p\ < r < s 
, r<-(s-p) 

2nr , otherwise 

defined for r > 0. 

The calculation is a little more intricate for the Stokes 
parameters. The light from each point on the star's disk 
must be linearly polarized in the tangential direction. Us- 
ing the angle x defined in Fig. [l], the Stokes parameters 
must therefore be F\j(r, x) = —P{ r ) sin 2% and Fq{t, x) — 
—P{r) cos 2x, for some function P(r) which we wish to re- 
cover (note that we use the unnormalised Stokes parame- 
ters, since the normalised ones have contributions to the 
noise from the intensity as well as the polarization mea- 
surements). Defining 

,.2ji 



r / cos2xA(r,x;s,c/))dx, 



(19) 



5. Polarization data inversion 

One strength of the Backus-Gilbert method is that we can 
do a lot of the analysis before we have any data. This gives 
us an understanding of the limitations of our analysis, 
and allows us to pick an optimal value for the smoothing 
parameter A. 

5.1. Theoretical limits on polarization inversions 

The Backus-Gilbert method is one of many inverse prob- 
lem methods in which one functional of the recovered so- 
lution, A, is minimised subject to regulation by another 
functional B. In this case, functional A is the width of 
the averaging kernel, and this is regulated by the vari- 
ance of the estimate. In other inverse problem methods, 
the functional A is some measure of the goodness of fit 
between the data and the forward problem, such as a x 2 , 
regulated by the demand that the solution be smooth, or 
that some non-linear functional of the solution, such as its 
negentropy, be minimised. 

The Backus-Gilbert scheme produces a one- 
dimensional family of solutions, parametrised by the 
smoothing parameter A. We can represent this as a solu- 
tion curve on a graph of standard deviation against of the 
recovered u against the chosen measure of the width of 
the resolution function: the resolution is inversely related 
to the width, so such a curve illustrates the trade-off 
between accuracy and bias in the recovered solution. 

In general, one usually considers intermediate values 
of A, around the turning point between the two extremes. 
Precisely which A one chooses depends on the relative im- 
portance of stability and resolution in the particular prob- 
lem. 

In the present case, we are interested in measuring the 
polarization at the limb, and it is clear from Fig. 2 that 
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0.05 0.1 0.15 0.2 0.25 

standard deviation 

Fig. 2. An example of a trade-off curve for the Backus-Gilbert inversion scheme. 



the resolution we need can only be obtained at the cost of 
a high standard deviation in the recovered value. If we try 
to increase the accuracy we end up no longer measuring 
the limb polarization proper, but rather a "blurred" value 
of the polarisation, smoothed by convolution with the av- 
eraging kernel. If the polarization is a maximum at the 
limb, the bias introduced by this smoothing will reduce 
the estimated polarisation: the more sharply the polariza- 
tion falls away from the limb, the greater this biasing will 
be. 

Looking at the low-A (ie, no-noise) limit, we find that 
the A(r, r ) function is sharply peaked at r = 1, so that in 
principle we can extract a well-resolved limb-polarization. 
In fact, the quality of this peak degrades substantially for 
r < 1: the method as presented cannot reasonably resolve 
polarizations on the disk. If we wished to recover these, 
we might use our knowledge of the kernels and include in 
the sum in Eqn. (fj]) only those kernels which are non-zero 
in [0, r]. Such a procedure would give no improvement for 
r = 1. 

As we increase A, the peak broadens, but the variance 
of the recovered u\(r) decreases. In Fig. [|we display the 
width A and variance B we obtain when we recover P(l) 



using the kernel Aq, for various values of A from 10 -2 
to 10 2 . Here we can clearly see the trade-off between the 
well-resolved but unstable recovery in the bottom right 
corner, and the stable but badly-resolved recovery in the 
top left. 

This diagram in a sense completes the Backus-Gilbert 
analysis of the inverse problem, as represented by the ker- 
nel in Eqn. (|2l]) , and we are now in a position to move on 
to invert real or simulated data. Before we can do that, 
however, we must decide what value of A to use. To make 
that decision, we must consider the level and approximate 
functional form of the polarization P(r), and use this to set 
the scale for the resolution and standard deviation we need 
to achieve. In turn, this fixes the number of data points 
n we require in our data and the value of the parameter 
A we must choose in our inversion. Despite the fact that 
we are invoking a particular model at this point in our 
analysis, we emphasise that this introduces no practical 
model dependence. We are using an approximate model 
purely to help us understand what counts as 'sufficiently 
stable' or 'sufficiently well resolved', and after this under- 
standing is gained the numbers we recover remain model 
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independent measurements, as opposed to any method of 
parameter fitting. 

Firstly, Chandrasekhar suggests that the limb polar- 
ization is of the order of P(l) = 0.1; we therefore need 
a variance which is at least as small as this, requiring 
log 10 Var[u] < —1. Secondly, if we are not to have an 
overly biased result, our resolution function must be nar- 
row compared with the width of the underlying func- 
tion P(r). The resolution we need is therefore of order 
Width[P(r)], with the width functional used in Eqn. (||). 
Taking P(r) ~ exp lOr as representative, we find we need 
a resolution better than log 10 Width [A] = —1.6. Compar- 
ison with Fig. ^ indicates that A = 1 and n 60 should 
therefore give us a satisfactory recovery of P(l). The max- 
imum accuracy achievable for a given resolution (or vice 
versa) can be read off from the graph. 

Fig. [5] shows the recovery of the limb polarization from 
three sets of simulated noisy data, as a function of the 
smoothing parameter. It is clear that the solution becomes 
more stable as greater smoothing is imposed. However, 
there comes a point where further smoothing does nothing 



to improve the solution, but merely degrades the resolu- 
tion of the recovery. 

5.2. The Algol system: a worked example 

As a more concrete illustration, consider the Algol system. 



Given the system parameters (from Sodcrhjelm 1980) we 
can analyse the observational prospects as in the previ- 
ous section. The Algol system is more complex than the 
spherically symmetric cases we consider, but our simpli- 
fied analysis will provide an upper bound on the resolution 
achievable with the real system. 



The graphs presented by Kemp et al. (1982) have 25 
data points in the eclipse phase. The same paper states 
that, in the Algol primary star, 50% of the polarized flux 
comes from an annulus on the limb of width less than 
0.005 R disc . 

Our Backus-Gilbert analysis of a spherically symmet- 
ric model system with the Algol orbital and luminosity 
parameters indicates that this limb polarization profile 
cannot be resolved in the Algol system . For the eclipse cov- 
erage reported in Kemp et al. ( 1983|) the minimum width 
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Fig. 4. The width of A(r,r'), log 10 A, versus log 10 B = log 10 Var[u], for a variety of values of A, and a variety of 
sample sizes n. Along each line, log 10 A rises from -2 in the bottom right corner, to +2 in the top left, with integer 
values of log 10 A plotted with a diamond on the two n — 60 lines. The resolution improves as we move down the plot, 
and becomes more stable as we move to the left; both the stability and the resolution improve as we use more data 
point s. T hese curves are for n = 20, 60, 112, as indicated; for the discussion of the line marked 'n = 60 (uneven)', see 
Sect. 5.3. The box in the bottom left encloses those combinations of A and n which will produce an estimator u which 



is adequately stable and well-resolved. 



-4min of the averaging kernel at the limb is 0.0199 i?disc, 
even in the zero-noise limit. 

This indicates that the limb polarization cannot be 
truly resolved with this data if the polarization profile is 
really as sharp as the stellar models cited by Kemp et al. 
predict. Thus, while the data does indicate the detection 
of limb polarization, it cannot reasonably be used (as was 
attempted by Wilson and Liou ( 1993 )) to make a quanti- 
tative estimate of the polarization profile. 

The precise dependence of the resolution on the num- 
ber of data points is beyond the scope of the present paper, 
but we note here that increasing the number of points from 
25 to as many as 1000 does not reduce the kernel width 
sufficiently (from 0.019 i?di sc to 0.0144 Ra\ sc )- For practi- 
cal purposes, then, the inadequate resolution is intrinsic 



to the Algol system and will not be alleviated by improved 
observations. 



5.3. Observational strategies 

Given a set of measurements fi at positions Sj, the Ba- 
ckus-Gilbert method can give us a well-controlled recov- 
ery of the underlying function u(r). We can do better than 
this, however, as the method can suggest how to improve 
our observational strategy to improve the resolution of the 
recovery. Clearly, increasing the total number of measure- 
ments we take, or improving the noise on each measure- 
ment, will improve the quality of our recovery. Equally 
clearly, simply binning the data (so that a ~ y/n) demon- 
strably improves nothing (this can be thought of as con- 
servation of information). 
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Fig. 5. The recovered values of P(l) for three realisations of noisy simulated data with the same parameters. The 
horizontal line is the correct value of P(l). 



Even if we assume, however, that we cannot change 
the number or quality of our data points (telescope time 
and equipment are limited, after all), we might still have 
the freedom to adjust the positions Sj (ie, the times) at 
which me make our measurements. 

In Eqn. (^), we are averaging the data points fi with a 
weight vector qi. Where qi is relatively large, therefore, /, 
will be smeared over several qi, or over a range of Sj, smear- 
ing out features in the kernel K(r; Sj). If we can identify 
prominent features in the vector q^, and cluster our mea- 
surements round the Sj they correspond to, we should be 
able to decrease the spread of A(r, r') and so increase the 
resolution of the recovery. 

Our simulations show that this rather informal argu- 
ment is valid for our case. The line in Fig. ^ captioned 
l n = 60 (uneven)' has the same number of data points as 
the line 'n = 60', but with the values Sj chosen so that the 
'data rate' in a band s = 1.85±0.15 is double that outside 
the band. This does not significantly improve the variance 
of the result (we are not gathering any more information 
than before), but it does noticeably improve its resolution. 



This improvement in our procedure corresponds to 
concentrating our measurements on the point around the 
beginning of the eclipse. We might have guessed that this 
would be a reasonable strategy to adopt, but the Backus- 
Gilbert method has justified our guess, and would substi- 
tute for a lack of intuition in a more obscure situation. 

6. Conclusions 

Our studies of simulated data show the fundamental limi- 
tations on the determination of limb polarization in eclips- 
ing binary stars. In particular, Fig. ^ indicates that limb 
polarizations of order a few percent are only just above 
the threshold of detectability, even in perfectly spherically 
symmetric, non-interacting binary systems. The situation 
will be worse in more complex systems. 

It is important to appreciate that the Backus-Gilbert 
method does not strictly estimate the polarization at a 
point on the stellar disc, but rather the polarization con- 
volved with the resolution function. To relate the results 
of the inversion to a particular model, it is necessary to 
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calculate the theoretical value of this convolution, which 
should be consistent with the A we have chosen and with 
all higher values of A - these represent coarser averages 
over the stellar disc. 

The bottom line is that one must take care in drawing 
conclusions about limb polarization from studies of eclips- 
ing binaries. It is clearly not possible to distinguish be- 
tween stellar atmosphere models on this basis if their pre- 
dicted limb polarizations differ by less than the maximum 
accuracy achievable. On the other hand, an appreciation 
of the issues raised in this paper will allow a meaning- 
ful determination of limb polarization, with reliable error 
estimates. 

Formally, the eclipsing binary problem is very similar 
to the gravitational microlensing problem. Indeed, part 
of the initial motivation for this work sprang from stud- 
ies of the microlens ing of extended sources. A future pa- 
per (Coleman et al. in preparation) will apply the inverse 
problem approach outlined in this paper to the use of mi- 
crolensing as a probe of stellar atmospheres. 

Acknowledgements. We would like to thank Dr. Richard K. 
Barrett for enlightening discussions of the Backus-Gilbert me- 
thod. IJC was supported by a PPARC studentship. 

Appendix A: General derivation of the 
Backus-Gilbert index 

In section 3, we described the Backus- Gilbert inverse in 
rather practical terms, taking care to relate the description 
to the quantities obtained in, and the concerns relevant to, 
real observations. Such a physical understanding of the 
method is essential if it is to be used properly, but we can 
obtain other insights into the method by reexamining it 
in a more formal wayQ. 

Eqn. (^) above describes an operator K : P — > D 
mapping an object from a source space P into a data 
space D. Including noise n S D, we have 



/ = Ku + n 



(A.l) 



for / G D and u G P. Here P is a real Hilbert space, 
parametrised by r, with a symmetric inner product 

(a\b) P = [ a{r)b(r)dr ,Va,6eP, 
Jo 

and D is a finite-dimensional Euclidean space with 
(a\b) D = y]aibj. 

i 

We wish to make an estimate u r S R of a single compo- 
nent of the object u, based on the data /. To this end, we 
wish to find a q G D (depending on r), such that 



(A.2) 

We thank the referee, Dr A Lannes, for suggesting this 
approach. 



We find this q as the solution of a minimisation problem. 
Introducing the adjoint operator K* : D — ► P, and as- 
suming E ({q\n)D) — 0, we have 

E (u r ) = (q\Ku) D = (K*q\u) P . (A.3) 

Recalling that (K*q\u)p = f (K*q)(r')u(r')dr', we see 
that (K*q)(r') G P can be identified with the averaging 
kernel A(r, r'), and that Eqn. (A.3) will be a good estimate 
of u r when q is transformed by K* into the basis vector 
e r G P corresponding to the component r of u. That is, 
Eqn. (A.3) would be exact if K*q = e r . The object K*q 
will instead be a linear combination of basis vectors 'close' 
to e r , and we can measure its 'scatter' around e r with the 
operator Q : P — > P such that Qx r > = (r — r') 2 x r > , Vav € 
P. Define 

A = (K*q\QK*q) P = (q\KQK*q) D = £ (ftWyfc, (A.4) 

defining the (self-adjoint) operator W = KQK* : D — > D. 
We may also define a measure of the stability of u r , 

B=(q\Sq) D , (A.5) 



by analogy with Eqn. (|11|), where the operator S G D 
is such that (et\Sej)D = S%j, where Sij is the positive 
definite noise covariance matrix. The demand that A(r, r') 
have unit area translates into the constraint (K*q\ l)p = I , 
where 1 G P is the all-1 vector in P. Writing R = Kl G D, 
this is equivalent to the constraint 



(q\R)D = 1, 



(A.6) 



restricting q to a hypersurface in D, with normal R. 

If we now introduce the functional c : D — > R, such 
that 

C(q) = \ \A{q) + \B(q)} = \ (q\ (W + XS)q) D , (A.7) 

the minimisation problem becomes that of finding the q 
which minimises c(q), subject to (q\R)jj — 1. Considering 
small variations e G D in the hyperplane (that is {e : 
(e\R) D = 0}, and defining the gradient Vc(q) = (W+XS)q 
(W + XS is self-adjoint since both W and S are), we have 



c(q + e) = c(q) + (e|Vc(g)) 



D 



XS)e] 



D ■ 



This is extremised at q\ such that (e| Vc(qa))d = 0, and 
is a minimum if (eKW^ + XS)e)pi > 0. The operator S 
is positive-definite by definition, and the operator W is 
positive-definite, since Q is. The operator c(q) is therefore 
minimised when (e|(W + XS)q\)n — 0,Ve, or 

(W + XS)q x = aR, 

for any a G R. Imposing the constraint Eqn. (|/ 
thus find 

(W + XSr'R 

qx {R\{W + XS)~V R D' [ 1 

from which we can obtain Eqn. (|l5|), from the definition 
of (-\-) D . 



(A.8) 

I), we 
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